Pipe vibration attenuation through internal damping and optimal design of vibro-impact systems

Pipelines periodically supported by rack structures (PPRs) are common in chemical and petrochemical plants, among others, and conventional tools such as dampers and hysteretic absorbers are commonly used to mitigate large vibrations in these systems. In this study, we explore two alternative strategies: (i) enhancing the attenuation rate of PPR vibrations through structural internal damping, and (ii) using nonlinear vibro-impact systems (VIS) to reduce seismic vibrations in a PPR. To shed light on the first strategy, we develop analytical dispersion relations for a PPR and show how damping can improve the mitigation capabilities of the periodic system. As for the second strategy, we consider a 9-node beam, i.e., a single span (SS) of a PPR equipped with a VIS, and combine the central composite design (CCD) and Kriging metamodelling to maximize dissipation energy and minimize the number of impacts. This multi-objective optimization problem aims to find the most effective design solution for the VIS in terms of gap and coefficient of restitution (COR). Additionally, we consider the stochastic nature of seismic input and the possible chaotic behavior of the VIS. To account for the sensitive variability of the number of impacts in seismic records, we perform incremental dynamic analyses and calculate fragility functions for various engineering demand parameters, including the number of impacts. We define a 3D surface for selecting the optimal gap-COR pair. When impacts occur, transient results can be chaotic, and we compute the largest Lyapunov exponents of a few representative trajectories.


Methods
Internal and external damping models in periodic PPRs. To improve the attenuation properties of the PPRs depicted in Fig. 2 both damping and periodic properties are employed. The supports of the pipelines depend on the flexural stiffness k v of the pipe rack's columns. In Fig. 2a, flexible springs support the periodic pipeline. Figure 2b depicts the model of a pipeline supported by rigid supports. Hereinafter, the models in Fig. 2a, b are called PPR f and PPR r , respectively. The damping models proposed to enhance vibration mitigation are linear and have different sources: internal material damping due to localized plastic deformation within the apparent elastic limit, see Fig. 2c; and the external viscous damping that is proportional to the forcing frequency, see Fig. 2d. Since the length dz → 0 , the damping is continuously distributed along the pipe. Therefore, internal and external damping forces can be represented by two different equations for the Euler-Bernoulli beam.
For the PPR f and the internal damping model the Euler-Bernoulli equation reads, where ρ, E, A and J are, respectively, the material density, the Young modulus, the area and the inertia of the beam cross section. The term w(z,t) is the beam transversal displacement. C i is the internal damping coefficient. Instead, for the externally damped PPR f , the Euler-Bernoulli equation reads, where C e is the external damping coefficient. In (1) and (2) the damping coefficients read, respectively, where ζ i/e is the modal damping ratio; i/e specifies internal or external damping and ω n defines the modal frequency. Note that the two coefficients C i/e have different physical meanings and are dimensionally different. C i in (3) comes from the assumption that structural internal damping does not entail plastic deformations in the cross section of the pipe 19 . In this regard, Kimball et al. 20 showed that, for metals subjected to cyclic stress, internal friction entails strains that remain below the elastic limit. C e in (4), instead, represents the common proportionality constant of a viscous damping model. The mass and the stiffness of the pipe are constant along its length. Since the systems are periodic, we consider (1) and (2) and apply the Floquet-Bloch theorem to the jth support, C e = 2ζ e,n ρA ω n ,  where ω defines the circular frequency. Full derivation of (5) is provided in the Supplementary information: Appendix A. About the PPR r , a dispersion relation between μ and ω is available in literature 15 and reads, Two-dimensional FEMs of the PPR f and PPR r have been modelled via Ansys APDL 19.0 with Euler-Bernoulli beam elements. A time-harmonic rotation ϕ i/p e i2πft was imposed as input at the left-end of the pipe, and the steady state response ϕ o/p (f) was read as the output rotation at the right-end. The system response in the frequency domain is evaluated as follows, where the FRF is the frequency response function expressed in dB. To foster the dynamics of a finite periodic system, the FEM consists of a 40 spans beam.
Optimization of the stochastic dynamic performance of the VIS. To analyse the performance of the VIS, nonlinear transient analyses were carried out for the equations of motion of a single unit cell of the PPR f , i.e., a SS beam supported by flexible springs. The aforementioned FEM was adopted to obtain the consistent mass matrix M and the symmetric stiffness matrix K, whereas the damping matrix C was computed with the proportional Rayleigh damping model. The discretized model of the SS equipped with a vibro-impact device is depicted in Fig. 3.
The interface of the SS has a two-sided displacement constraint, such that the node g-5 impacts two rigid bumpers. In this case, we directly model the VIS and let the system of equations of motion to remain fully linear. The coefficient of restitution (COR) is defined according to Newton's law 26 : where v r is the relative velocity of two colliding bodies, and + and − mean post-impact and pre-impact, respectively. Since the bumpers are fixed to the ground, relative velocities coincide with absolute velocities. Note that Eq. (10) does not need further nonlinearities to be included. However, to take into account the typical limitations of the aforementioned instantaneous model 26,34-36 , very small time steps of the order of 10 -5 have been used; thus, the simulations showed a limited dependence from time steps 36 .
The system of equations of motion read, where M, C and K are respectively the mass, damping and stiffness matrices, r is the influence vector and u g is the ground displacement; the dot represents the derivative wrt the time. When impact occurs, the g node is (5) ψ i/e cosh 2 (µ i/e )+χ i/e cosh(µ i/e )+η i/e = 0, Figure 3. Single span (SS) discretized into 9 nodes that vibrate vertically and rotate; f indicates the nodes that do not experience the impact, whereas g-5 has a two-sided displacement constraint. where R g is the impact force vector. Equation (12) is divided in two different equations, one for the non-impacting nodes f and one for the impacting node g. Since the impact force is unknown, we consider the first equation of the system in (12), as To obtain the structural response at the ith step of the algorithm, Eq. (12) is employed when the pipe and the bumper are not in contact. Conversely, if the contact occurs, (13) is used. A flowchart that describes the algorithm implemented in MATLAB for the nonlinear impact model is shown in Fig. 4.
At a certain instant t d the energy dissipated by linear viscous damping is computed as follows: whilst the energy dissipated by the impacts reads, where n d is the number of impacts occurred up to t d , and ẋ − and ẋ + represent the velocity vector before and after impact, respectively. Once the seismic event is extinguished, free decay oscillation occurs in the structure www.nature.com/scientificreports/ up to a certain time t r at which the rest condition is fully restored, and the number of impacts is n im . The total input seismic energy reads: Therefore, at time instant t r , all the input seismic energy has been dissipated by means of viscous damping and impacts. In these conditions, the design parameters COR and gap that optimize the seismic performance of the VIS are sought by a multi-objective optimization problem based on two objective functions defined as follows, where n im is the number of impacts, E imp and E tot are defined by (15) and (16), respectively. The purpose of installing bumpers next to the pipe is to dissipate the largest amount of energy through impacts whilst limiting the number of impacts, to prevent damage to the pipe. Therefore, the objective function O 1 must be maximized (equivalently, − O 1 must be minimized), while O 2 is to be minimized. To reduce the number of transient analyses, we combine CCD and Kriging modelling with the explanatory variables COR, gap, and Sa (T 1 ), and response variables O 1 and O 2 . More precisely, Sa (T 1 ), i.e., the intensity measure (IM), represents the spectral acceleration at the first period T 1 of the system. The 15 yellow points of Fig. 5 indicate the computer-experimental data set-sampling points-defined by the CCD. The ranges of values of the explanatory variables are shown in Table 1.
Note that, to take into account the stochastic nature of the seismic input, the IM Sa(T 1 ) has been included in the CCD. The system is nonlinear and sometimes chaotic, see Supplementary Information-Appendix B; therefore, a certain variability of the system response wrt the input is expected. As indicated in Table 1, Sa (T 1 ) is bounded by the mean minus/plus the standard deviation of the selected records. Then, to evaluate O 1 and O 2 for the 15 points selected within the CCD, seismic transient analyses are carried out. Each ground motion, see Table 2, is scaled at the relevant value of Sa(T 1 ) and thus, O 1 and O 2 are calculated as the mean of the values obtained by employing all the seismic records described hereinafter. After running all the 15 × 12 analyses needed for CCD, a Gaussian polynomial regression between the 15 mean values by means of the Kriging metamodeling provided by UqLab 43 has been performed. It has been assumed that the model output is a realization of a Gaussian process defined as the joint distribution of the prediction and the true model response 43 .
Then, the Pareto front is used to provide the optimal values of the multi-objective optimization problem by seeking the nondominated solutions among the O 1 and O 2 quantities of (17) and (18) evaluated through the Kriging metamodel. To each nondominated solution of the Pareto front will correspond an optimal triplet (Sa(T 1 ), gap, COR). Therefore, the optimal values of gap and COR are connected to a certain value of Sa(T 1 ); nonetheless Sa(T 1 ) is not a true design parameter since it characterizes a (seismic) stochastic process.   Fig. 3, a set of twelve natural records with a 2% probability of exceedance in 50 years, are employed; that is, relevant to safe shutdown events (SSE). The selection of the natural seismic records follows the principle sketched in Fig. 6. Both the mean spectrum and the mean spectrum plus one standard deviation of the selected accelerograms match with the uniform hazard spectrum (UHS) of a specific site, Priolo Gargallo, Sicily in Italy, in a least-square sense. More precisely, let us consider s 0 the target spectrum value vector, that is, the UHS, and evaluate as the spectra matrix of the n a accelerograms. One can define a vector of n a selection coefficients, α, where each element can only take a binary value of 1 or 0, and the sum of the elements is equal to n s , i.e., the number of accelerograms to be selected. Thus, the vector α that satisfies is sought. The selection is performed with all possible combinations of the n s accelerograms among a set of n a records. This operation allows to preserve full seismic hazard consistency of the site and minimize the recordto-record variability yet considering the dispersion of the records about the mean spectrum. Table 2 reports the set of accelerograms and their main characteristics.
Finally, we perform the fragility assessment of the SS. Along this line, the fragility function F DM (IM) is defined as the probability that the node g-5 of Fig. 3 reaches or exceeds some damage measure (DM) for a given ground motion with IM = im. In particular, the DM is connected to the number of impacts n im , i.e., the engineering demand parameters (EDP), that exceeds a certain threshold. Typically, F DM (IM) is assumed to follow a lognormal distribution, and hence reads,  (20) is calculated for all the feasible n im = EDP corresponding to the Pareto front. In fact, the Pareto front highlights the optimal design parameters gap and COR. Therefore, each optimal value of the Pareto front corresponds to a fragility surface (F DM (IM), IM, EDP). The corresponding volume V j under the jth surface is computed as follows where j indicates the jth optimal solution of the Pareto front. According to the law of total probability, (23) is proportional to the total probability of exceedance of all the edps. Then, one should perform the fragility assessment for all the j = {1,…,n} Pareto front solutions and select the optimal couple gap-COR that minimizes (23). In this work, an application of this procedure is shown for three optimal couples. When the pillars' stiffness is considerably higher than the pipe's flexural stiffness, the pillars act as rigid supports wrt the flexural behaviour of the pipe. Thus, the PPR r is modelled as a beam supported by simple supports. The dispersion diagrams of the internally damped PPR r , are calculated via Eq. (6) and are depicted in Fig. 7. Figure 7 shows the dispersion diagrams (μ−f) for seven values of material damping ζ i of Eq. (3). The yellow bands in Fig. 7a Fig. 7 reveal common properties of linear periodic systems. The blue curve represents the undamped pipe and clearly defines the bandgaps in the regions where Re(μ) ≠ 0 and, consequently, Im(μ) = 0 or Im(μ) = π, unveiling the existence of pure evanescent waves. Instead, for damped pipes with nonzero values of ζ, Re(μ) ≠ 0 in the overall frequency domain, since attenuating oscillatory waves can always be observed for all frequencies. However, the damping is included in the main structure, and clearly enlarges the attenuation rate Re(μ). This effect is more evident as the frequency range increases. Figure 7c shows the FRF of (9), evaluated with the FE software. A good agreement is observed between Fig. 7c and the dispersion curves of Fig. 7a, b, in terms of passbands, bandgaps and attenuation rates.

Results
Let us now consider the PPR f depicted in Fig. 2b. The pipe rack is modelled by a spring-mass oscillator that matches the first lateral mode of the rack. Iqbal et al. 17 demonstrated the effectiveness of such an approximation. The relevant dispersion curves are calculated via (5) and plotted in Fig. 8. Fig. 8a, b show the dispersion relations of the PPR f , where a local resonance is observed around the first natural frequency of the pipe rack ω pr,1 = 4.45 Hz. The first bandgap, that initially ranged from 0 to 5.25 Hz, is now divided into two different bandgaps. Generally, bandgaps can be tuned with local resonances, whereas the Bragg scattering induces bandgaps that are constrained by the periodic system dimensions. Nevertheless, in this paper the PPR f 's first natural frequency belongs to the first bandgap frequency range, and therefore no additional bandgaps appear. Conversely, Fig. 8a, b depict a zoom on a narrow passband that opened due to the pipe rack's resonance. Finally, one can notice a good agreement between the numerical FRF of Fig. 8c and the analytical dispersion diagrams of Fig. 8a, b. Again, the damping induces an attenuation rate that increases with the frequency.
The combination of Eq. (7) with (8) allows for plotting the dispersion diagrams and the bandgaps relevant to the external damping case. Also in that situation, favourable bandgap zones can be obtained; for brevity, they are not shown or commented on. (20) F DM (IM) ≡ P EDP > edp|IM = im , www.nature.com/scientificreports/ Optimization results. As anticipated in a previous subsection, the Pareto front S1 depicted in Fig. 9 pinpoints the nondominated solutions among all the possible realizations evaluated with the Kriging model. This Pareto front is called S1 to avoid confusions with Fig. 12.
The surrogate Kriging metamodel generates O 1 and O 2 as 3D-arrays for the parameters gap, COR and Sa(T 1 ). Therefore, the Pareto front highlights the optimal triplets (Sa(T 1 ), gap, COR), and the optimal values of gap and COR corresponds to a certain optimal value of Sa(T 1 ). Nonetheless, Sa(T 1 ) is characterized by randomness and influences the optimal solutions. For clarity, the three surfaces that depict the values of O 1 wrt three optimal Sa(T 1 ), equal to 0.6 g, 1.3 g and 2 g, are depicted in Fig. 10.
Note that the objective function O 1 defined in (17), must be maximized. Nonetheless, we plot -O 1 in Fig. 10 and then seek its minimum values. The surfaces O 2 of (18) refer to the optimal Sa( 1 ) equal to 0.6 g, 1.3 g and 2 g in Fig. 11.
The black circles in both Figs. 10 and 11 represent the Pareto front of Fig. 9.
To underline the sensitivity of the optimization problem to the seismic input Sa(T 1 ), Fig. 12 shows the Pareto fronts corresponding to several levels of Sa(T 1 ).
It is clear that the optimal solutions differ significantly as the severity level of the seismic input varies. For clarity, Fig. 13 depicts all the Pareto fronts of Fig. 12 along with the relevant surfaces O 1 and O 2 .
One can clearly note that the optimization is sensitive to the seismic input Sa(T 1 ); therefore, the inclusion of the IM in the CCD is justified. A few more general considerations arise from the optimization results. Both the objective functions O 1 and O 2 depend more on the COR than on the gap. In particular, optimal values of COR are found in the whole range 0-1, whereas large gap values are optimal for strong earthquakes, say Sa(T 1 ) equal to 2 g, 2.7 g, and 3.4 g. For lower values of Sa(T 1 ), the optimal gaps cover the whole range 40-80 mm.
Fragility assessment and selection of the optimal solution. To perform the fragility assessment, three optimal solutions among those of the Pareto front in Fig. 9 are randomly taken; and the IDA are calculated for the eight ground motions indicated with the * in Table 2. Figure 14 reports the results of the IDA and the relevant 2D fragility functions.
The impact phenomenon may generate chaotic motion, as described in the Supplementary Information-Appendix B; hence the IDAs result in responses that are not unique for a certain EDP = n imp . Vamvatsikos et al. 52 described the occurrence of multiple capacity points, and recommended to handle this ambiguity by an ad hoc,     Finally, and for the sake of completeness, the reader can appreciate the speed of the algorithm sketched in Fig. 4 by using Optimal solution #1. In this respect, Table 3 reports the CPU time required to simulate the response of the optimal VIS to the records listed in Table 2 using MATLAB. The simulations run on a processor with a clock speed of 3.00 GHz and 4 cores.
One can observe that the CPU time employed to solve the system of Eq. (11) is effective in relation to the event duration.

Discussion and outlooks
In this paper, the enhanced attenuation properties of two periodic damped pipelines coupled with pipe racks (PPRs) have been shown; and the relevant results have been confirmed by the FE software Ansys on a 40-spans PPR. Two damping models were proposed in view of vibration mitigation, i.e., internal material damping, and external viscous damping. A local resonance was observed around the first natural frequency of the PPR ω pr,1 = 4.45 Hz. The damping clearly enlarges the attenuation rate Re(μ), and this effect is more evident as the frequency range increases. As a result of the analyses, we have found a good agreement between the frequency response function (FRF) of the FE models and the analytical band structures.
Then, a generic single span (SS) of the PPR equipped with a vibro-impact device was considered; and due to the nonlinearities, a design optimization procedure was conceived and carried out. The procedure aims to maximize the dissipation energy and to minimize the number of impacts. Both objective functions resulted to be more sensitive to the COR than to the gap; the response surface was generated by the Kriging metamodel with a Gaussian 2nd order-polynomial regression. It was found that the metamodels were endowed with sharp curvatures wrt to COR. Optimal values of COR were found in the whole range 0-1; large gap values appear to be optimal for strong earthquakes, say Sa(T 1 ) equal to 2 g, 2.7 g, and 3.4 g, whereas optimal gaps are found in the whole range 40-80 mm for lower values of Sa(T 1 ). Moreover, we have found that the optimal solution with gap = 64 mm and COR = 0.49 is the one that minimizes the probability of exceeding the damage states (DM) connected to the engineering demand parameter (EDP), that is, the number of impacts. The incremental dynamic analysis (IDA) resulted in curves that are not unique for a certain EDP, i.e., noisy, mainly due to the chaotic www.nature.com/scientificreports/ behaviour of the VIS. In addition, the IDA curves displayed significant record-to-record variability: for example, in Fig. 14a, the EDPn im = 3000 was reached for values of Sa(T 1 ) equal to about 1 g, 1.6 g, 1.7 g, 2.4 g, 2.5 g, 5 g, 5.3 g and 7 g. This result, however, must be imputed to the frequency content of the records and the nonlinear response of the vibro-impact that exhibits a strong dependency upon frequency and amplitude of excitation, and sometimes chaos. To investigate this last issue, we have excited the single span (SS) with periodic loadings for a restricted range of frequencies in the vicinity of ω pr,1 , and have reported the relevant bifurcation diagrams in the Supplementary Information -Appendix B. When impact does not occur, the bifurcation diagram shows that the system is linear and periodic. Conversely, the impact activates higher modes of vibrations, and nonperiodic solutions that can be found in the bifurcation diagram. Therefore, three of these trajectories have been investigated, and the largest Lyapunov exponents were found to be higher than zero, thus indicating the presence of divergence and chaos.
In conclusion, we have shown how to evaluate the safest design solution for a nonlinear dissipation system despite the (seismic) stochastic nature of the loading. Nevertheless, some challenges remain undisclosed. The defined dispersion curves concern the linearly damped system without impact-induced dissipation. Nonlinear waves, in fact, distort as they propagate and change their original shape along the periodic medium; therefore, numerical transient analyses of a multiple-span pipe rack with vibro-impact devices may reveal unintended periodic features of the impacting repetitive structure.